/*
 * BeadPos.h
 *
 *  Created on: Jan 26, 2011
 *      Author: snirgaz
 */

#ifndef BeadPos_H_
#define BeadPos_H_
#include "HelpLibs/def.h"

using namespace std;

template<typename E>
// A CRTP base class for Vecs with a size and indexing:
class VecExpression {
public:
	typedef array<double, SIM_DIM> container_type;
	typedef container_type::size_type size_type;
	typedef container_type::value_type value_type;
	typedef container_type::reference reference;

	size_type size() const {
		return static_cast<E const&>(*this).size();
	}
	value_type operator[](size_type i) const {
		return static_cast<E const&>(*this)[i];
	}

	E& get_ref() {
		return static_cast<E&>(*this);
	}
	E const&get_ref() const {
		return static_cast<const E&>(*this);
	}
};

// The actual Vec class:
class ParticlePos: public VecExpression<ParticlePos> {
	container_type _data;
public:
	reference operator[](size_type i) {
		return _data[i];
	}
	value_type operator[](size_type i) const {
		return _data[i];
	}
	size_type size() const {
		return _data.size();
	}
	ParticlePos() {
	}
	bool inBoxLimits() {
		bool flag = true;
		for (int d = 0; d < SIM_DIM; d++){
			double pos=_data[d];
			if (( pos > 0.5) || (pos < -0.5)) {
				flag = false;
				cout << pos << endl;
			}
		}
		return flag;
	}
	void rotate() {
		// double rot_param[2][2] = { 1, INV_TAN_T, 0, 1 / sin(PI / 3) };
		// Rotate
		_data[0] = _data[0] + _data[1] * INV_TAN_T;
		_data[1] = _data[1] * INV_SIN_T;
	}
	void rotate_inv() {
		//double rot_param[2][2] = { 1, cos(PI / 3), 0, sin(PI / 3) };
		// Rotate
		_data[0] = _data[0] + _data[1] * COS_T;
		_data[1] = _data[1] * SIN_T;
	}
	void cyclic() {
		for (int d = 0; d < this->size(); d++) {
			if (_data[d] > .5)
				_data[d] -= 1;
			if (_data[d] < -.5)
				_data[d] += 1;
		}
	}
	double sumOfSqures() const {
		double sum = 0;
		for (int d = 0; d < SIM_DIM; d++)
			sum += pow(_data[d], 2);
	#ifdef ROT
		sum += data_[0] * data_[1];
	#endif
		return sum;
	}
	double norm2() const {
		return sqrt(this->sumOfSqures());
	}
	// Construct from any VecExpression:
	template<typename E>
	ParticlePos(VecExpression<E> const& vec) {
		E const& v = vec.get_ref();
		for (size_type i = 0; i != v.size(); ++i) {
			if (v[i] > .5)
				_data[i] = v[i] - 1;
			else if (v[i] < -.5)
				_data[i] = v[i] + 1;
			else
				_data[i] = v[i];
		}
	}
	template<typename E>
	void operator+=(VecExpression<E> const& vec) {
		E const& v = vec.get_ref();
		for (size_type i = 0; i != v.size(); ++i) {
			if (v[i] > .5)
				_data[i] += v[i] - 1;
			else if (v[i] < -.5)
				_data[i] += v[i] + 1;
			else
				_data[i] += v[i];
		}
	}
	template<typename E>
		void operator-=(VecExpression<E> const& vec) {
			E const& v = vec.get_ref();
			for (size_type i = 0; i != v.size(); ++i) {
				if (v[i] > .5)
					_data[i] -= v[i] - 1;
				else if (v[i] < -.5)
					_data[i] -= v[i] + 1;
				else
					_data[i] -= v[i];
			}
		}
};
// The actual Vec class:
class Force: public VecExpression<Force> {
	container_type _data;
public:
	reference operator[](size_type i) {
		return _data[i];
	}
	value_type operator[](size_type i) const {
		return _data[i];
	}
	size_type size() const {
		return _data.size();
	}
	Force() {
	}
	// Construct from any VecExpression:
	template<typename E>
	Force(VecExpression<E> const& vec) {
		E const& v = vec.get_ref();
		for (size_type i = 0; i != v.size(); ++i) {
			_data[i] = v[i];
		}
	}
	void fill(value_type d) {
			for (size_type i = 0; i != size(); ++i) {
				_data[i] = d;
			}
		}
	double sumOfSqures() const {
		double sum = 0;
		for (int d = 0; d < size(); d++)
			sum += pow(_data[d], 2);
#ifdef ROT
		sum += data_[0] * data_[1];
#endif
		return sum;
	}
	template<typename E>
			void operator+=(VecExpression<E> const& vec) {
				E const& v = vec.get_ref();
				for (size_type i = 0; i != v.size(); ++i) {
						_data[i] += v[i];
				}
			}
	template<typename E>
			void operator-=(VecExpression<E> const& vec) {
				E const& v = vec.get_ref();
				for (size_type i = 0; i != v.size(); ++i) {
						_data[i] -= v[i];
				}
			}
};
template<typename E1, typename E2>
class VecDifference: public VecExpression<VecDifference<E1, E2> > {
	E1 const& _u;
	E2 const& _v;
public:
	typedef ParticlePos::size_type size_type;
	typedef ParticlePos::value_type value_type;
	VecDifference(VecExpression<E1> const& u, VecExpression<E2> const& v) :
			_u(u.get_ref()), _v(v.get_ref()) {
		// assert(u.size() == v.size());
	}
	size_type size() const {
		return _v.size();
	}
	value_type operator[](size_type i) const {
		return _u[i] - _v[i];
	}
};

template<typename E1, typename E2>
class VecSum: public VecExpression<VecSum<E1, E2> > {
	E1 const& _u;
	E2 const& _v;
public:
	typedef ParticlePos::size_type size_type;
	typedef ParticlePos::value_type value_type;
	VecSum(VecExpression<E1> const& u, VecExpression<E2> const& v) :
			_u(u.get_ref()), _v(v.get_ref()) {
		//   assert(u.size() == v.size());
	}
	size_type size() const {
		return _v.size();
	}
	value_type operator[](size_type i) const {
		return _u[i] + _v[i];
	}
};

template<typename E>
class VecScaled: public VecExpression<VecScaled<E> > {
	double _alpha;
	E const& _v;
public:
	typedef ParticlePos::size_type size_type;
	typedef ParticlePos::value_type value_type;
	VecScaled(double alpha, VecExpression<E> const& v) :
			_alpha(alpha), _v(v.get_ref()) {
	}
	size_type size() const {
		return _v.size();
	}
	value_type operator[](size_type i) const {
		return _alpha * _v[i];
	}
};

// Now we can overload operators:

template<typename E1, typename E2>
VecDifference<E1, E2> const operator-(VecExpression<E1> const& u,
		VecExpression<E2> const& v) {
	return VecDifference<E1, E2>(u, v);
}

template<typename E1, typename E2>
VecSum<E1, E2> const operator+(VecExpression<E1> const& u,
		VecExpression<E2> const& v) {
	return VecSum<E1, E2>(u, v);
}
template<typename E>
VecScaled<E> const operator*(double alpha, VecExpression<E> const& v) {
	return VecScaled<E>(alpha, v);
}

///////////////////////////////////////////////////////
struct Bead {
	// world line number
	int n_;
	// time bead slice
	int RWSlice_;
	// time bead total path
	int WLSlice_;
	// Bead Position;
	ParticlePos beadPos_;
	// Jump Bead Position;
	ParticlePos jumpPos_;
	//
	bool originalJump_, proposedJump_;
	// Is Odee
	bool isOdd;
	void set(int n, int RWSlice, int WLSlice) {
		n_ = n;
		RWSlice_ = RWSlice;
		WLSlice_ = WLSlice;
	}
	static bool beadsOrder(Bead const &first, Bead const &second) {
		return (first.RWSlice_ < second.RWSlice_);
	}
	bool isJump() {
		return proposedJump_;
	}
	bool isOriginalJump() {
		return originalJump_;
	}
};
typedef boost::shared_ptr<Bead> BeadP;
typedef std::vector<Bead> BeadVec;
// Beads Functions
bool beads_comp(Bead first, Bead second);
bool oddBead(Bead b);

#endif /* BeadPos_H_ */
